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Abstract 

In  this  paper  we  investigate  oscillations  of  a  dynamical  system  containing  passive  dynamics 
driven  by  a  positive  feedback  and  how  spatial  characteristics  (i.e.  symmetry)  affect  the  ampli¬ 
tude  and  stability  of  its  nominal  limit  cycling  response.  The  physical  motivation  of  this  problem 
is  thermoacoustic  dynamics  in  a  gas  turbine  combustor.  The  spatial  domain  is  periodic  (pas¬ 
sive  annular  acoustics)  which  are  driven  by  heat  released  from  a  combustion  process,  and  with 
sufficient  driving  through  this  nonlinear  feedback  a  limit  cycle  is  produced  which  is  exhibited 
by  a  traveling  acoustic  wave  around  this  annulus.  We  show  that  this  response  can  be  controlled 
passively  by  spatial  perturbation  in  the  symmetry  of  acoustic  parameters.  We  find  the  critical 
parameter  values  that  affect  this  oscillation,  study  the  bifurcation  properties,  and  subsequently 
use  harmonic  balance  and  temporal  averaging  to  characterize  periodic  solutions  and  their  sta¬ 
bility.  In  all  of  these  cases,  we  carry  a  parameter  associated  with  the  spatial  symmetry  of  the 
acoustics  and  investigate  how  this  symmetry  affects  the  system  response.  The  contribution  of 
this  paper  is  a  unique  analysis  of  a  particular  physical  phenomena,  as  well  as  illustrating  the 
equivalence  of  different  nonlinear  analysis  tools  for  this  analysis. 


1  Introduction 

Thermoacoustic  instabilities  in  gas  turbines  develop  when  acoustics  in  a  combustor  couple  with  an 
unsteady  heat-release  in  a  positive  feedback  loop.  Thermoacoustic  modeling  and  control  is  well- 
studied  for  axially  extended  combustion  chambers,  as  in  [1,  2,  3,  4],  where  the  acoustic  to  heat- 
release  coupling  is  dominated  by  longitudinal  acoustic  modes.  Because  different  instability  regimes 
occur  at  different  operating  conditions,  recent  attention  has  focussed  on  thermoacoustic  modeling 
in  combustion  chambers  with  annular,  or  cylindrical  geometries  [5,  6,  7,  8,  9]. 

Control  of  thermoacoustic  oscillations  is  a  rich  field  due  to  the  complexity  of  the  dynamics  as  well 
as  the  prominence  of  both  land  and  air-based  jet  engines.  These  high-energy  devices  operate  in  a  wide 
range  of  operating  conditions,  all  of  which  are  highly  nonlinear  due  to  turbulence,  combustion,  and 
other  extreme  conditions.  The  oscillations  lead  to  compromised  performance,  high  noise  levels,  or 
catastrophic  engine  damage.  Current  control  means  (see  [10])  include  avoiding  operating  conditions 
which  exhibit  large  oscillations,  additional  dissipation  including  acoustic  dampers  or  resonators  [11], 
or  in  some  cases  active  feedback  control  [12,  13,  14].  The  first  option  is  detrimental  in  terms  of 
marketing/engine  use  while  redesign  using  dampers  and  resonators  takes  time  and  adds  weight. 
Active  control  is  challenging  due  to  actuator  limitations  including  harsh  conditions,  high  temporal 
frequencies,  and  the  limited  amount  of  control  authority  a  finite  number  of  actuators  offers. 

With  these  limitations  in  mind,  any  approach  that  opens  the  operating  envelope,  and  does  not 
add  weight  or  significant  complexity  and  redesign  is  a  valuable  solution.  In  this  paper  we  study 
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such  a  solution;  introducing  precise  spatial  variations  (asymmetry)  in  a  specific  mean  property  of 
the  dynamics  which  directly  affects  the  amplitude  of  limit  cycle  oscillations. 

Approaches  to  analyze  thermoacoustic  oscillations  are  abundant,  ranging  from  complex  CFD 
simulation  to  reduced  order  analysis.  The  reduced  order  analysis  often  includes  spatial  approxima¬ 
tion  of  the  first  principle  dynamics  resulting  in  a  system  of  ordinary  differential  equations.  This 
spatial  approximation  can  be  performed  many  ways  including  an  element  wise  approach  or  by  a 
modal  decomposition.  If  a  modal  approach  is  used,  the  model  is  often  of  very  low  order  (since  the 
thermoacoustic  instability  is  typically  tied  predominantly  to  one  or  two  acoustic  modes)  while  the 
element  approach  is  used,  the  ODE  system  will  result  in  numerous  coupled  oscillators.  We  mention 
this  here  because  although  our  approximation  in  this  study  is  modal,  resulting  in  two  coupled  modes 
(oscillators),  the  tools  used  here  align  with  those  used  in  to  study  synchronized  aspects  of  multiple 
coupled  oscillators. 

Interaction  of  nonlinear  oscillator  populations  has  been  a  rich  topic  for  hundreds  of  years.  From 
pendulum  clock  synchronization,  to  fire-ffy  hashing,  to  the  interaction  of  spiking  neural  activity,  the 
study  of  the  behavior  of  these  systems  continues  to  offer  insight  into  either  amplification  or  control 
of  synchronized  behavior  in  interconnected  complex  systems.  In  [15]  coupled  oscillators  in  a  ring 
using  uniform  distribution  of  natural  frequencies  were  analyzed  seeking  stable  conditions  including 
synchrony  and  oscillator  death.  Our  work  seeks  similar  information  while  the  natural  frequencies 
(acoustic  wavespeed)  of  the  oscillator  population  are  perturbed  by  a  spatial  pattern  instead  of  a 
statistical  distribution.  Similarly  traveling  wave  solutions  for  a  ring  of  nonlinear  oscillators  was 
investigated  in  [16]  while  a  completely  different  form  of  network  connection  was  investigated  in  [17]. 

This  study  differs  from  typical  synchronization  studies  predominately  because  we  investigate  a 
modal  response  of  the  oscillators  instead  of  seeking  conditions  on  n  individual  oscillators.  In  addi¬ 
tion,  most  synchronization  studies  investigate  the  behavior  of  oscillators  which  are  self  excited  in 
their  uncoupled  state  while  here,  the  coupling  itself  motivates  a  limit  cycle.  Finally,  most  often 
synchronization  studies  investigate  binary  dynamics  that  exhibit  either  on-off  or  integrate-and-fire 
responses  while  with  our  system  different  amplitude  solutions  are  prevalent  (phase  only  representa¬ 
tions  do  not  suffice).  Aside  from  these  differences,  the  methods  used  here  align  significantly  with 
those  used  in  the  aforementioned  efforts. 

The  organization  of  this  paper  is  as  follows;  we  begin  by  developing  a  reduced  model  of  the 
dynamics  starting  with  first  principles.  The  equilibrium  of  this  model  is  then  investigated  using 
numerical  bifurcation  analysis  tools.  Following  this,  the  periodic  equilibriuma  are  determined  using 
two  methods;  Averaging,  and  Harmonic  Balance.  The  results  of  these  two  methods  are  compared 
and  the  limit  cycle  properties  are  investigated  using  both  approaches  illustrating  that  altering  the 
symmetry  in  a  specific  parameter  always  reduces  the  limit  cycle  amplitude  and  eventually  stabilizes 
the  system.  It  is  also  shown  that  the  oscillations  during  this  imposed  asymmetry  remain  stable. 

2  Modelling 

In  this  section  we  briefly  describe  the  thermoacoustic  model  used  in  the  current  analysis.  In  an  annu¬ 
lar  combustion  system,  flow  passes  down  the  length  of  an  annulus,  and  is  eventually  mixed  with  fuel. 
When  this  mixture  reaches  a  series  of  flames  distributed  around  this  annulus,  it  reacts/combusts. 
This  combustion  process  both  reacts  to  the  acoustics  (i.e  how  the  acoustics  bring  the  fuel  mixture 
to  the  flame),  while  at  the  same  time  it  drives  the  acoustics  with  its  heat  release.  Figure  1,  presents 
a  schematic  of  the  annular  combustor  and  feedback  interconnection  between  the  acoustics  and  heat 
release. 

The  fundamental  structure  and  justification  of  all  assumptions  of  the  transport  equations  for  this 
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Figure  1: 


Schematic  of  the  Annular  Domain  with  the  First  Pair  of  Fourier  Modes 


system  is  available  in  the  works  of  Culick  (i.e.  [18]).  The  transport  equations  are: 

0  (1) 

-Vp  (2) 

-pV  ■  u  +  q.  (3) 

where  p,  u,  and  e  are  density,  velocity,  and  internal  energy  per  unit  volume,  p,  and  q  are  pressure  and 
a  volumetric  heat  release  contribution  and  the  spatial  Laplacians  describe  cylindrical  coordinates. 

Using  a  series  of  assumptions  relating  to  the  boundary  conditions  which  can  be  found  in  [22] 
results  in  a  pair  of  partial  differential  equations  on  the  annular  domain: 


|  +  v.M  = 

du  ^ 

+  pu-Vu  = 
at 

de 

p^  +  pu-Ve  = 


dp 


due 

dt 

due 

de 


=  -Cp  +  q, 


(4) 

(5) 


where  0  is  the  spatial  rotational  coordinate,  ^  is  a  damping  constant,  a  is  the  acoustic  wavespeed, 
and  9  is  a  driving  heat  release  mechanism. 


2.1  Heat  Release 

The  equations  (4,  5)  are  self  contained  with  the  exception  of  the  heat  release  function.  The  form  of  the 
heat  release  function  (g)  is  the  driving  component  in  the  feedback  loop.  Because  of  the  complexity  of 
combustion  dynamics,  harsh  conditions,  high  noise  levels,  and  lack  of  adequate  sensing  apparatus,  an 
analytical  form  is  not  yet  available.  The  heat  released  from  combustion  is  most  frequently  modeled 
as  a  function  of  acoustic  velocity  that  contains  a  saturation  like  quality  as  in  [1,  3,  19,  20,  21].  With 
this  in  mind  we  denote  the  heat  release  contribution  as: 

q  =  a[u]+a[u^]  (6) 

where  cr  is  a  linear  destabilizing  effect,  and  a  is  a  parameter  that  introduces  the  limiting  effect  (a  <  0 
is  chosen). 

2.2  Asymmetry  Parameter 

Prior  to  discretizing  the  equations  into  a  modal  representation,  the  acoustic  wavespeed  is  parame¬ 
terized  to  have  a  spatial  preference.  A  perturbation  is  added  to  the  wavespeed  in  the  form  of  the 
second  spatial  harmonic.  The  reason  for  this  choice  is  that  using  linear  analysis  in  a  previous  study 
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[22],  it  was  found  that  altering  the  wave-speed  in  this  pattern  effected  the  stability  the  most.  The 
resulting  spatially  perturbed  acoustic  wavespeed  is: 

a  =  ao  +  cos(20)  (7) 

where  oq  is  a  nominal  wave-speed  and  S  is  the  asymmetry  parameter  with  second  harmonic  prefer¬ 
ence. 

2.3  Discretization 

The  discretized  equations  are  obtained  by  projecting  the  rotational  coordinate  onto  spatial  Fourier 
modes.  The  first  two  spatial  modes  are  retained;  sin(0), cos(0).  A  spatio-temporal  response  in 
these  coordinates  is: 

p{Q,t)  =  {psin{t)  Sm{Q)  +  Pcos{t)  cos{Q))  (8) 

u{Q,t)  =  {us^n{t)  sin{Q)  +  Ucos{t)  cos{Q))  (9) 

where  Psin,Pcos  and  itsm,  Ucos  are  the  time  dependent  modal  amplitudes  which  become  the  states  of 
our  dynamical  system.  A  schematic  of  these  modes  on  the  annular  domain  is  presented  in  Figure  2. 


Figure  2:  Schematic  Representation  of  the  Spectral  Modes  of  the  Annular  Combustor 


The  relative  amplitude  and  phase  of  the  temporal  coefficients  impact  the  presence  of  either 
standing  or  traveling  waves.  A  standing  wave  is  characterized  as  a  amplitude  that  depends  on 
position  (i.e.  a  response  that  has  clear  nodes  and  antinodes),  while  a  traveling  wave  retains  a 
constant  amplitude  with  space.  These  wave  motions  have  reciprocal  relations  in  that  traveling 
waves  can  be  described  using  combination  of  standing  waves  and  vice  versa  [23] .  We  are  particularly 
interested  in  traveling  wave  solutions  and  keep  in  mind  that  using  sin(0)  and  cos(0)  as  a  basis  for 
the  dynamics,  a  traveling  wave  exists  when  the  amplitudes  of  the  temporal  coefficients  are  equal 
with  a  phase  difference  of  |  (and  solutions  that  are  trigonometrically  symmetric  to  this). 

When  using  the  first  two  Fourier  modes  as  a  basis,  the  cubic  nonlinearity  acts  through  the  velocity 
component  and  is  projected  onto  the  modes  as  follows: 

o;  /3  /tt  \ 

Qszn  =  ^  J  sin(0)  (Usinit)  sin(0)  -h  u^osit)  cos(0))^  de  =  a  f  j  (10) 

(x  /3  \ 

dcos  =  J  cos(0)  (Usinit)  sin(0)  -F  Ucosit)  cos(0))^  dO  =  a  f  {uliri  +  uLs)  j  (H) 

and  projecting  the  assumed  modal  solution  onto  (4,  5)  results  in  the  system  of  ODE’s: 
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To  get  these  equations  in  the  form  coupled  oscillators,  we  take  [usin,  Mcos]  as  [xi,X2\  which  results 
in  the  following  second  order  coupled  nonlinear  differential  equations: 


x'l  +  Cx[  +  {al-5)xi  =  X2{al-6)  (^a  -  a  {xj  +  xj)^^  (12) 

X2  +  (^2  +  {al  +  S)x2  =  -xi{al  +  S)  ^tr  +  a  (a;?  +  (13) 

From  this  representation,  the  architecture  of  the  system  becomes  clear.  Once  projected  onto  the 
two  modes,  the  dynamics  are  two  stable  oscillators  with  skew  symmetric  coupling  due  to  the  heat 
release. 

For  the  remainder  of  this  study,  we  will  investigate  the  behavior  of  system  (12,  13)  specifically 
focusing  on  two  parameters;  a  which  is  the  skew  symmetric  coupling  parameter  from  heat  release, 
and  5  which  is  the  asymmetry  parameter.  We  begin  by  studying  behavior  of  the  global  equilibrium 
and  continue  by  investigating  the  periodic  equilibria  and  their  stability.  For  the  analysis  to  follow, 
we  will  use  the  values  in  Table  1  as  a  set  of  nominal  parameters. 


Table  1:  Nominal  Parameters 


Oq 

c 

a 

a 

6 

1.0 

0.2 

0.5 

-0.0018 

0.0 

3  Equilibrium  Analysis 

The  system  (12,  13)  contains  one  equilibrium  point  at  the  origin,  and  in  this  section  we  investigate  the 
stability  of  this  equilibrium  with  respect  to  both  the  linear  coupling  parameter,  and  the  asymmetry 
parameter.  When  the  asymmetry  parameter  (5)  and  coupling  term  (cr)  are  zero,  evaluation  of  the 

Jacobian  at  the  origin  reveals  a  set  of  co-located  stable  eigenvalue  pairs  at  ^  ±  a/— Joq  +  . 

When  the  coupling  is  increased  from  zero  these  complex  eigenvalues  become: 

Ai,2  —  - 

A3.4  = 

It  is  evident  that  with  sufficiently  large  values  of  positive  or  negative  coupling  the  eigenvalues  break 
apart,  moving  left  and  right  in  the  complex  plane.  Physically,  the  case  with  cr  =  0  represents  zero 
combustion  (the  eigenvalues  are  purely  acoustic),  and  the  case  with  nonzero  a  portrays  the  case 
with  the  driving  combustion  process.  With  sufficient  heat  addition,  one  pair  of  complex  eigenvalues 
eventually  crosses  the  imaginary  axis.  The  value  of  the  coupling  parameter  such  that  the  eigenvalues 
cross  the  imaginary  axis  is: 

^  crit  —  (^^) 

Oo 

From  this  we  see  that  instability  occurs  when  the  coupling  from  positive  feedback  exceeds  the 
acoustic  damping  normalized  by  the  nominal  wavespeed.  When  the  coupling  from  heat  release  is 
further  increased  the  eigenvalues  continue  to  become  more  unstable,  the  oscillations  grow  and  are 
eventually  limited  by  nonlinearity.  Figure  3  shows  results  from  the  numerical  bifurcation  analysis  tool 
AUTO  illustrating  the  Hopf  Bifurcation  leading  to  instability  without  any  modification  of  symmetry 
(the  system  parameter  is  oscillation  amplitude). 

The  stabilizing  effect  of  adding  asymmetry  to  the  combustion  dynamics  is  performed  by  investi¬ 
gating  the  behavior  of  the  system  under  variations  of  the  asymmetry  parameter  <5.  Taking  a  nonzero 
positive  value  for  cr  which  insures  oscillations  in  the  dynamics,  we  perform  eigenvalue  analysis  in  a 


C+^J  C^-4ag±4-yUT^ 
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Bifucation  Diagram 


Figure  3:  Bifurcation  Solutions  Representing  the  Destabilizing  Effect  of  the  Parameter  a 


similar  way.  The  results  show  that  increasing  asymmetry  provides  a  stabilizing  mechanism,  and  the 
critical  value  that  returns  the  eigenvalues  to  the  stable  half  of  the  complex  plane  is: 


QoVQpg^  -  C 

\/l  + 


(16) 


Similar  numerical  bifurcation  results  of  the  system  with  its  symmetry  altered  showing  its  stabi¬ 
lizing  effect  (Figure  4). 


Bifucation  Diagram 


Figure  4:  Bifurcation  Solutions  Representing  the  Stabilizing  Effect  of  the  Parameter  6 

The  focus  of  this  study  is  to  generate  both  analytical  and  numerical  results  which  are  similar 
to  the  schematics  of  Figures  3  and  4.  Specifically,  we  investigate  the  amplitude  and  stability  of  the 
limit  cycle  when  cr  >  Gcrit  and  0  <  (5  <  Scrit- 
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4  Periodic  Solutions 


To  analyze  the  bifurcation  branches  and  periodic  behavior  of  the  combustion  system,  we  rely  on 
the  characteristic  that  the  dynamics  are  nearly  linear.  The  linearity  of  the  system  is  addressed  by 
rearranging  the  terms  of  the  nonlinear  system  (12,  13)  into  a  standard  oscillator  form  and  entering 
in  the  nominal  parameters.  The  equations  become: 

x"  +  xi  =  -0.2a;;  +  6xi  +  0.5(-(5  +  l)a:2  +  0.0024(5  -  l)a;;a;2  +  0.0024(5  -  l)xl  ,  . 

x"  +  X2  =  -0.2a;'2  -  5a;2  +  0.5(-5  -  l)a;i  +  0.0024(5  +  l)xixl  +  0.0024(5  +  l)a;f 

With  small  variation  in  the  asymmetry  parameter  5  the  system  remains  closely  linear  in  its  behavior. 
This  is  also  confirmed  by  the  mostly  sinusoidal  response  of  equations  from  time  integration  (i.e. 
simulation  with  an  ODE  solver).  Another  point  of  interest  which  is  evident  in  (17)  is  that  the  terms 
on  the  right  hand  side  consist  of  driving  from  heat  addition  and  dissipation  from  damping  in  the 
acoustics.  It  is  a  balance  of  these  two  mechanisms  that  is  the  manifestation  of  the  periodic  response 
of  the  system. 

Since  the  dynamics  of  the  system  appear  to  first  order  as  those  of  a  linear  oscillator  with  small 
perturbation,  we  have  a  few  options  for  approximate  analysis  of  the  dynamics.  Below  we  present 
two  approaches  averaging  (in  two  different  coordinate  systems),  and  harmonic  balance  to  obtain  the 
slow  flow  dynamics  of  the  trajectories. 

4.1  Averaging 

The  first  approach  to  obtain  the  slow  flow  of  equations  (12,  13)  includes  the  use  of  time  averaging 
(see  [24]).  We  perform  this  approximation  using  two  different  basis,  polar  coordinates,  and  using 
action  angle  variables.  We  show  that  the  averaged  results  are  identical  using  either  approach. 

The  averaging  approach  includes  making  a  harmonic  assumption  of  the  limit  cycle  response 
exploiting  a  time  scale  assumption  on  the  assumed  solution  (i.e.  variation  of  parameters),  and 
integration  in  time  (averaging).  What  results  is  the  slow  flow  (ie  amplitude  and  phase  ODE’s).  We 
we  will  investigate  the  equilibria  of  this  reduced  system  both  with  and  without  assymetry. 

4.1.1  Averaging  in  Polar  Coordinates 

We  begin  by  making  the  assumption  that  the  response  is  harmonic  with  slowly  varying  amplitude 
and  phase  parameters: 


Xi{t)  =  Ri{t)Cos{ujt  +  (j)i{t))  (18) 

where  i  relates  to  the  two  oscillators  describing  the  amplitudes  of  the  sine  and  cosine  acoustic 
modes.  For  this  study,  we  are  concerned  with  solutions  where  the  period  of  each  oscillation  is  equal 
and  stationary.  Differentiation  in  time  results  in: 

x'{t)  =  —ujRi{t)Sin{ujt  +  4'i{t))  —  4)^(1) Ri{t) Sin{ujt  +  (j)i{t))  +  i?((t)  cos{ujt  +  (j)i{t))  (19) 

where  (.)'  represents  differentiation  with  respect  to  time.  The  assumption  that  the  amplitude  and 
phase  of  the  assumed  solution  vary  slowly  with  time  reveals  that  their  time  derivative  has  negligible 
contribution  to  the  velocity  states  which  is  equivalent  to  solution  of  the  ODE  using  variation  of 
parameters[26].  With  this  in  mind,  we  assume  : 

0  =  -(j)'i{t)Ri{t)  sm{Lut  +  (flit))  +  i?'(f)  cos(wt  +  (fiit))  (20) 

which  results  in  the  velocity  state  as: 

x'it)  =  —ujRiit)  siniut  +  (fiit))  (21) 

The  assumed  solutions  (18),  and  (21),  and  the  time  derivative  of  (21)  are  substituted  into  the 
equations  of  motion.  Using  this  equation  and  (20)  results  in  two  equations  and  two  unknowns  for  each 
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oscillator.  Isolating  the  time  derivatives  of  the  unknowns  results  in  an  ODE  describing 

the  slow  evolution  of  the  amplitude  and  phase  of  the  limit  cycle.  This  differential  equation  system 
will  contain  the  parameters  of  the  system  (cr,  a,  oq,  ■Co)  the  states  of  the  slow  flow  {Ri{t) ,  ,  and 

resonant  terms  (i.e.  cos(a;t), sin(a;t)).  The  next  step  is  to  remove  the  resonant  terms  by  averaging 
over  one  period  of  oscillation.  Before  doing  this  we  change  the  coordinates  to  consider  relative  phase 
variables: 


V'+  =  5  {0i  +  62) ,  Ip-  =01-92 
or  6»i  =  ’^  +  V’+,  02  =  1^+  - 


(22) 


and  perform  the  averaging: 


dipi 

dt 


T 


f{R,ip,t)dt 


(23) 

(24) 


where  T  =  ^,  and  the  function  f{R,ip,t)  contains  information  from  both  oscillators  due  to  the 
coupling,  and  the  overbar  denotes  averaged  quantities.  The  resulting  system  is: 


^  TT  ]  64aoHi^2  I  ]  ^ 

-f  32i5i?ii?2+cos('0_  )  +  — 

32aoi?ii?2 

16Cao^i+sin(i/j_  )  (aQ  — i5).R2  (16(7+3\/¥c^(^i+3.R|)  ) 

~  ~  32ao 

16Cao^2+sin(i/j_  )  (aQ+(5)i?i  (l6fT+3v^a(3i?i+i?2)) 

.  ^0^ 


(25) 


Studying  the  equilibrium  points  of  (25)  reveals  the  amplitude  and  phase  of  the  limit  cycle  when 
the  limit  cycle  in  the  system  exists  due  to  sufficient  heat  release  (i.e.  a  >  <Jcrit)-  We  will  study  these 
equilibria  and  how  they  are  effected  by  assymmetry  in  Section  5. 


4.1.2  Averaging  in  Action  Angle  Coordinates 

To  perform  the  averaging  in  action  angle  coordinates,  we  first  define  the  canonical  transformation: 

^i{t)  =  •\/2Ji(t)  cos  {0i{t)) ,  xi{t)'  =  ^/2Jl{t)s\a{0l{t))  (26) 

X2{t)  =  \/2J2{t)  cos  {02{t)) ,  X2{t)'  =  yj2J2{t)sm{92{t))  (27) 

where  again  as  designated,  the  action  and  angle  variables  are  functions  of  time.  The  second 

derivative  in  time  is  found  from  differentiation  of  the  above  relations  (i.e.  x'(  =  -^^)  and  af¬ 
ter  substituting  the  coordinate  transforms,  we  have  two  equations  and  four  differential  variables 
( Ji,  J2,  01, 6*2).  Two  additional  equations  are  realized  by  the  constraint  equations  (again  from  vari¬ 
ation  of  parameters): 


/ 

ij 


/ 

2 


(28) 


for  nature  of  the  assumed  quantities  in  26,  and  27. 

After  these  substitutions,  we  can  solve  for  the  differential  variables  of  the  system  ( J{,  J2,  0(,  0^) 
resulting  in  a  fourth  order  ordinary  differential  equation  system.  Again,  we  transform  the  phase 
variables  to  the  form  in  Equation  22  resulting  in  a  differential  equation  system  (J(,  J2,  V'+i '*/'-)■ 
this  coordinate  system,  the  phase  variable  ip+  acts  like  time  and  hence  we  average  the  dynamics 
over  this  variable: 


0 


fdp}+ 


(29) 


The  averaged  system  in  Action  Angle  coordinates  is  then: 


j  /  A(fio  +  l)  +cos(^ip_  ^  ^S(T Ji  + J2))  ((  A  ~  A)fio+‘5(-/i  +  J2)) 

^  32y^yA 

7/  _  e  cos(i/'-)(8(T+9V7rQ(ji+J2))((A+A)o§+<5(ji- J2)) 

^_-6+ 

J'l  =  |\/^  fsin  [ip-]  (flo  “  '^)  {8a-  +  S^ra  (Ji  +  3J2))  -  8C\/Ji^ 

^2  =  \\/ ^2  ^sin  {'tpj)  (oq  +  |5)  \/ji  (8cr  +  Sa/^o  (3Ji  +  ^2))  -  Spy/j^j 

4.2  Harmonic  Balance 

The  variable  coefficient  harmonic  balance  method  provides  another  approximation  to  the  limit  cycle 
dynamics  and  to  develop  it  we  begin  with  the  assumed  harmonic  solution  (18),  differentiate  as 
necessary  and  insert  into  the  equations  of  motion  (12,  13).  Similar  to  the  use  of  variation  of 
parameters,  the  coefficients  of  the  assumed  harmonic  solution  are  considered  slowly  varying.  Because 
of  this,  as  in  [25],  the  second  order  derivatives  of  the  amplitude  and  phase  coefficients  of  the  assumed 
solution  are  considered  negligible.  By  collecting  the  sine  and  cosine  terms  of  each  equation  (the 
coefficients  of  the  resonant  terms)  we  obtain  two  differential  equations  for  each  oscillator.  Similar 
to  the  slow  flow  system  derived  from  Averaging  this  ODE  describes  the  slow  evolution  of  amplitude 
and  phase  of  the  assumed  solution.  We  may  again  solve  for  the  time  derivatives  of  the  unknown 
coefficients,  and  the  equilibriums  will  describe  the  limit  cycle  conditions. 

If  we  consider  the  resonant  case,  the  period  of  oscillation  will  be  related  to  the  wavespeed  and 
therefore  the  natural  frequency  of  the  assumed  solution  will  equal  the  wave  speed  (w  =  ao).  After 
including  this  into  the  equations  of  motion,  we  have  the  slow  flow  system: 


R'i  = 
i?' = 

V''-  = 


f,'  = 


3v^a(3k2C+2klao)(oo-5)fl2fl?  +  16C(<5-2ao)fli  +  (k2C+2klao)(ao-5)fl2(9v^ai?2  +  16(T) 

16(C^+4a§) 

-3VSct(2klao-3k2C)(ao+(5)i?ifl2+16C(2ao+<5)i?2-(2klao-k2C)(ao+5)fli(9Wrai?i+16cr) 

16(c^+4af) 

9v^a(klC+2k2oo)(ao+5)fi)+2(3Vtra(6k2a^+kl5C)A2+8(T(klC+2k2oo)(ao+5))fli 

16(C2+4a2)flii?2  ' 

645aoi?2-Ri  +  (2k2ao  — kl^)(aQ  — (5)i?2 
16(C2+4a2)i7i_R2 

-9V7ra(klC+2k2ao)(ao+i5)fl)-2(3Wtaao(6k2<5+klCao)i?2+8cr(klC+2k2ao)(ao+5))fli 

32(C+ial)RiR2  ’  ’ 

32C^aoi?2Ai  +  (2k2ao-klC)(ao-<5)fl2(9v^aA2  +  16o') 

32(C2+4a2)i?ii72 


where  kl  =  sin 'ijj_  and  k2  =  cosip-. 


(31) 


5  Comparison  of  Methods  and  Numerical  Results 

In  this  section  we  compare  the  results  of  the  different  approximation  methods.  We  first  present 
the  dynamics  when  the  asymmetry  parameter  is  set  to  zero  (the  value  of  the  equilibria,  and  the 
associated  phase  space  of  the  reduced  system) .  A  comparison  is  also  presented  for  the  amplitudes  of 
the  limit  cycle  with  varying  asymmetry.  It  is  shown  that  the  three  different  approximations  to  the 
full  nonlinear  system  are  in  agreement  with  each  other  and  that  the  asymmetry  parameter  always 
reduces  the  amplitude  of  the  limit  cycle. 

5.1  Reduced  Phase  Space  Comparison 

When  the  asymmetry  parameter  is  set  to  zero  the  dynamics  are  symmetric  and  therefore  the  am¬ 
plitude  of  the  two  acoustic  modes  are  identical.  This  offers  a  means  to  investigate  the  system  on  a 
two  dimensional  plane. 
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For  the  case  with  averaging  in  polar  coordinates,  we  let  Ri  =  R2  =  R  and  5  =  0  in  Equation 
(25)  and  obtain  the  reduced  system: 


=  0 

V’-  =  cos  {'4’-)  (iSy/naR'^  +  16cr)  qq  (32) 

H(l2v^a^^  +  16a-)  sin(i/)_  )aQ  +  16CHao 
^  ~  32a;; 

The  equilibrium  points  of  32  are  found  by  setting  time  derivatives  to  zero.  The  nontrivial  real 
solutions  are  presented  in  Table  2,  with  the  numerical  substitution  using  the  nominal  parameters. 

In  the  approximation  using  averaging  in  action  angle  coordinates  we  investigate  the  nominal 
symmetric  case  by  letting  Ji  =  J2  =  J,  <5  =  0  in  the  differential  equation  system  (30).  The  resulting 
equations  are: 


44+ =h  {-al-  1) 

4j'_  =  1  cos  {4j-)  (8cr  +  ISy^aJ)  ag  (33) 

J'  =  I VJ  ^■\/J  (8a  +  12y/7raJ)  sin  ('0-)  Oq  —  8CVjJ 

The  equilibrium  points  of  33  are  presented  in  Table  2.  Note  that  the  conversion  of  i?  « 
provides  exact  agreement  in  these  results  between  averaging  in  polar  and  action  angle  coordinates 
(and  the  harmonic  balance  for  that  matter). 

For  the  case  where  Harmonic  Balance  is  utilized  as  an  approximation,  we  investigate  the  sym¬ 
metric  case  by  letting  Ri  —  R2  =  R  and  (5  =  0  resulting  in: 


R'{t) 

R'{t) 

C  = 

0V  = 


cos('i/?_)C  — 16C+4(3v^ct-R^+4(T^  sin('0_)ao) 
8(C^-l-4ag) 

RaQ(^4:(3^/naR^ sin(T/i_)ao  — C((9\/^ck-R^+8o-)  cos('?/j_  )  +  16)) 
=  8(C=+4ag) 

(gV^afl^+So-)  cos(t/i_  )aQ 
2(C"+4a2) 

Cao(4C— (3\/w«fi^-l-4a')  sin(»/>_)ao) 

4(C2+4ag) 


(34) 


In  this  system,  the  term  0+  is  a  decoupled  non-equilibrium  state  that  acts  similar  to  time. 
Therefore  we  omit  this  differential  equation  and  perform  the  steady  state  analysis.  The  nontrivial 
real  equilibria  of  the  symmetric  system  are  presented  in  Table  2. 


Table  2:  Equilibria  for  the  Symmetric  Case  (Comparison  of  Different  Methods) 


- 

Averaging  (Polar) 

Averaging  (Action- Angle) 

Harmonic  Balance 

Amplitude 

D  2-v/±C  — (TOO 

—  1  , _ 

7r4  ^y3o'ao 

=  {17.10,11.19} 

j _ |_  2(CTf^<3,Q  j 

=  {62.68,146.27} 

D  2-v/±^— crao 

^  —  1  , _ 

TT  4  ^/3o^aQ 

=  {17.10,11.19} 

Phase 

'0-  =  ±f  ,±f 

0-  =  ±f  ,±f 

41-  =  ±f  ,±f 

As  seen  in  the  table  there  is  agreement  between  all  three  approximate  methods.  Additionally, 
since  the  system  has  been  reduced  to  a  second  order  system,  we  can  investigate  its  phase  space 
graphically.  The  phase  space  of  the  averaged  system  (32)  is  presented  in  Figure  5  (the  phase  space 
for  all  three  methods  is  indistinguishable  graphically,  so  for  brevity  we  present  the  vector  field  from 
averaging  in  polar  coordinates).  From  this,  we  find  that  the  equilibrium  point  R  =  11.19,  =  —  f 

is  stable,  while  all  other  equilibria  are  unstable,  which  was  confirmed  by  evaluation  of  the  Jacobian 
at  these  equilibria. 

5.2  Amplitude  Comparisons 

Now  that  we  are  familiar  with  the  approximate  system  when  the  asymmetry  parameter  is  zero,  we 
investigate  the  response  of  the  limit  cycle  amplitude  when  the  symmetry  of  the  acoustic  wavespeed 
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Figure  5:  Phase  Space  of  the  Nominally  Symmetric  System  Reduced  to  a  2D  Field 


is  modified  on  the  annular  domain.  In  this  section  we  provide  a  comparison  of  oscillator  amplitude 
(12,  13)  using  three  different  methods  including  averaging  in  polar  and  action  angle  coordinates,  as 
well  as  the  harmonic  balance  with  nonzero  asymmetry.  In  each  case,  fixed  points  of  the  ordinary 
differential  equations  were  found  by  time  integration  of  the  system  of  equations,  retaining  the  final 
point  when  the  dynamics  have  settled.  In  all  cases,  the  nominal  system  parameters  are  those 
described  in  Section  2,  while  the  symmetry  parameter  is  varied  as  the  independent  variable.  Figure 
6  illustrates  the  amplitude  of  the  limit  cycle  for  each  acoustic  mode  (each  oscillator  in  (12,  13)). 


Figure  6:  Mode  Amplitudes  as  a  function  of  the  Asymmetry  Parameter  (the  Action  Angle  Results 
are  scaled) 

In  this  figure,  the  influence  of  the  asymmetry  parameter  (S)  presents  an  attenuating  influence 
on  the  amplitudes  of  the  oscillations  ultimately  quenching  the  response.  In  addition,  we  notice  that 
these  amplitudes  decrease  at  different  rates  which  occurs  because  when  altering  the  symmetry  of  the 
system,  the  dynamics  prefer  one  spatial  mode  over  the  other,  and  the  corresponding  state  variable 
behaves  accordingly.  It  is  also  clear  that  the  solutions  agree  quiet  well  between  the  original  and 
approximate  dynamics. 

The  above  plots  trace  the  stable  solutions  of  the  system  in  both  its  full  and  approximate  state. 
We  are  also  interested  in  the  effect  of  asymmetry  on  the  unstable  solutions.  Below,  in  Figure  7 
trace  both  the  stable  and  unstable  branches  as  the  asymmetry  is  increased.  For  this  we  simply  use 
a  gradient  search  on  the  right  hand  sides  of  the  approximate  differential  equations  (seeking  zero) . 

This  figures  shows  that  when  5  =  0  and  =  i?2  there  exists  two  values  (11.19,17.10)  corre¬ 
sponding  to  the  stable  and  unstable  equilibria.  With  an  increase  in  5  the  stable  solution  begins  to 
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Limit  Cycle  Amplitude  vs.  Assymetry  Parameter 


Figure  7:  Mode  Amplitudes  as  a  function  of  the  Asymmetry  Parameter  for  Three  Different  Approx¬ 
imations  (note  all  approximations  are  nearly  equivalent,  and  the  Action  Angle  Results  are  scaled) 


decrease  in  amplitude  (  i?i  and  i?2  at  different  rates)  and  eventually  reaches  zero.  For  the  unstable 
solution  which  starts  at  17.10  increasing  6  also  effects  the  amplitudes  (at  different  rates)  but  it  is 
evident  that  these  solutions  do  not  interact  with  the  stable  branch. 


5.3  Orbit  Stability 

As  we  have  shown,  symmetry  effects  amplitudes  while  to  be  sufficiently  applicable  for  engineering 
application,  these  new  orbits  must  be  stable  to  account  for  unmodeled  characteristics  or  external 
influences.  Therefore,  it  is  necessary  to  assess  the  stability  of  the  periodic  solutions  of  this  system 
under  influence  of  the  asymmetry  to  characterize  the  robustness  of  the  use  symmetry  in  the  system 
to  reduce  oscillation  amplitudes.  To  perform  this  analysis,  the  slow  flow  equations  developed  from 
the  averaging  method  in  action  angle  coordinates  (30)  are  linearized  about  the  stable  equilibria  in 
Figure  7.  The  eigenvalues  of  these  fixed  points  is  presented  in  Figure  8. 


Eigenvalues  of  Periodic  Orbits  for  Various  Symmetry  Parameters 


Real 


Figure  8:  Eigenvalues  of  the  Periodic  Orbits  with  the  influence  of  Asymmetry 
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The  figure  above  illustrates  that  the  orbits  remain  stable  under  all  conditions  where  the  slow 
fiow  equations  are  valid.  The  behavior  of  these  dynamics  are  such  that  all  solutions  are  stable  and 
real  when  the  symmetry  parameter  is  set  to  zero.  Upon  increasing  this  parameter,  two  eigenvalues 
break  apart  and  become  complex  pairs.  Clearly  these  eigenvalues  are  associated  with  the  amplitude 
evolution  dynamics  (as  verified  by  investigating  the  eigenvectors),  while  the  third  eigenvalue  remains 
real  and  is  associated  with  the  differential  phase  between  the  two  oscillator  dynamics.  This  particular 
eigenvalue  approaches  instability  as  the  amplitude  equilibrium  approaches  zero  but  remains  stable. 


6  Summary 

In  this  paper  we  have  modeled  thermoacoustic  dynamics  in  an  annular  combustor  and  reduced 
these  dynamics  to  a  coupled  oscillator  system  including  a  parameter  for  the  spatial  symmetry  of  the 
passive  dynamics.  We  have  shown  that  the  nonlinear  coupling  due  to  heat  release,  always  imposes 
an  adverse  effect  on  the  system,  resulting  in  a  limit  cycle  exhibited  by  a  traveling  acoustic  wave.  We 
have  also  shown  that  spatial  perturbation  in  the  symmetry  of  the  acoustic  wavespeed  reduces  the 
amplitude  of  the  limit  cycle  (in  a  stable  way)  eventually  stabilizing  the  system.  These  results  have 
been  concluded  using  three  different  nonlinear  analysis  tools. 
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